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1.  Introduction 


A  large  amount  of  research  interest  has  focused  on  the  multiscale  problem  involving  atoms 
and  continua.  It  is  widely  accepted  that  many  effects  on  the  continuum  germinate  at  the 
atomic  level.  Events  such  as  fracture,  fatigue,  and  inelastic  material  response  can  be  traced 
back  to  phenomena  in  the  atomic  structure.  Moreover,  fabrication  of  nanoscale  devices  in 
mass  quantities  will  likely  come  from  nanopatterning  techniques  (e.g.,  see  [1]).  Creating 
patterns  at  the  nanometer  scale  will  likely  involve  unforeseen  effects  when  coupled  to  me¬ 
chanical  loads.  Therefore,  a  computational  mechanics  method  is  necessary  that  can  couple 
disparate  scales  -  one  scale  in  which  the  boundary  conditions  are  applied  and  the  other  in 
which  atoms  reside. 

Nanopatterning  has  recently  become  popular  for  creating  organized  nanostructured  materi¬ 
als,  which  are  those  that  may  be  defined  as  materials  whose  structural  elements,  clusters, 
crystallites  and  molecules  have  dimensions  in  the  1-  to  100-nm  range.  Experimental  methods 
can  be  designed  to  exploit  the  symmetries  and  repetitiveness  of  periodic  nanostructures  to 
control  the  environment  in  which  ultra-thin  films  and  surfaces  are  created.  Periodic  struc¬ 
tures  are  often  desireable  because,  when  atoms  and  electrons  are  confined  within  nanoscale 
semiconductors  and  metal  clusters,  unique  properties  arise.  Such  materials  are  formed  under 
a  different  paradigm  of  “bottom  up”  instead  of  conventional  “top  down”  techniques.  Figures 
1  and  2,  respectively,  illustrate  some  recent  work  of  experimental  and  numerical  nanopat¬ 
terning.  In  performing  rigorous  analysis  of  such  systems,  the  impact  of  mechanical  stresses 
and  strains  on  the  evolution  and  characteristics  of  the  atoms  must  be  modeled  carefully.  This 
requires  a  technique  for  translating  mechanical  information  at  the  boundaries,  e.g.,  clamps, 
pin  joints,  and  loads,  down  to  the  atomic  level.  Moreover,  to  ensure  that  the  mechanical 
nature  of  the  problem  stays  consistent  with  the  atomistic  problem,  the  information  contained 
at  the  atomic  level  must  be  transferred  back  up. 
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(a)  Array  of  nitrogen-  (b)  Cobalt  grown  on  the  (c)  Iron  deposited  on  the 

modified  copper  nanodots  copper-nitrogen  nanodot  sur-  copper-nitrogen  nanodot  sur- 

(dark  patches)  and  clean  face.  face, 

copper  lines  (bright  areas). 

Figure  1.  Scanning  tunneling  microscopy  (STM)  images  of  experimental 
nanopatterned  systems  [2], 

Methodologies  for  linking  a  continuum  to  an  atomistic  domain  can  be  found  in  literature  as 
early  1971  [4].  Finite  element  (FE)  methods  were  later  employed  in  Mullins  and  Dokain- 
ish  [5]  using  a  numerically  decoupled  domain  approach  with  spatially  overlapping  atomistic 
and  continuum  regions.  A  review  of  some  of  these  methods  can  be  found  by  Cleri  et  al.  [6]. 
Among  these  early  analytic  and  computational  studies,  frequent  issues  regarding  the  treat¬ 
ment  of  the  interface  arose,  which  were  primarily  handled  through  creative  use  of  kinematic 
constraints.  For  example,  Tadmor  et  al.  [7]  developed  an  FE-based  formulation,  the  so-called 
quasicontinuum  method.  Similar  efforts  were  made  through  the  so-called  handshaking  or 
coupling-of-length-scales  (CLS)  method  by  Broughton  et  al.  [8]  by  increasing  the  atomic 
resolution  to  account  for  electron  degrees  of  freedom  via  the  tight-binding  (TB)  method. 
The  dynamic  problem  was  studied  with  a  generalized  scaling  approach  in  coarse-grained 
molecular  dynamics  (CGMD)  by  Rudd  and  Broughton  [9]  to  better  handle  the  propagation 
of  waves  through  the  atomistic-FE  interface  and  the  FE  far  field. 


Figure  2.  Monte  Carlo  numerical  simulation  of  ion  erosion  of  platinum  (111) 
surface.  The  darker  diamond  region  is  the  periodic  cell  and  dark¬ 
ened  pit-regions  are  areas  of  greater  depth,  (a)-(d)  illustrate  re¬ 
moval  of  mono-layers  by  ion  erosion,  each  of  greater  depth  [3]. 

Multiscale  methods  such  as  these  have  traditionally  been  limited  mainly  to  localized  regions 
of  interest.  For  example,  the  applications  to  which  these  methods  have  been  applied  involve 
small  sets  of  dislocations  and  cracks  and  limited  analyses  of  their  mutual  interactions.  The 
localized  regions  on  which  these  simulations  are  run  typically  span,  at  most,  several  microns 
because  of  the  bottleneck  imposed  by  a  direct  interface  between  the  continuum  region  and 
atomistic  region.  To  ensure  compatibility,  kinematic  constraints  are  used  to  tie  together 
the  equations  and  disparate  length  scales  across  this  interface.  Driving  the  resolution  of  the 
discretized  continuum  down  to  the  atom  scale  to  accommodate  the  interface  intrinsically 
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restricts  the  size  of  the  continuum  and  leads  to  smaller  overall  dimensions  of  the  problem. 
This  can  only  be  overcome  by  larger  use  of  computer  resources  when  dealing  with  problems 
with  larger  dimensions. 

The  asymptotic  expansion  homogenization  method  has  been  widely  studied  by  applied  math¬ 
ematicians  for  many  years.  Numerous  texts  on  the  basic  theory  can  be  found  in  the  liter¬ 
ature,  for  instance,  by  Bensoussan  et  al.  [10],  Sanchez-Palencia  [11],  and  Bakhvalov  and 
Panasenko  [12].  Yet,  despite  the  prolific  research  in  the  field,  no  attempts  have  been  docu¬ 
mented  for  extending  the  technique  to  atoms. 

Homogenization  may  be  particularly  suited  for  nanopatterned  systems  because  of  the  use 
of  periodicity  and  asymptotics  in  the  assumptions,  perhaps  even  more  so  than  traditional 
materials  that  employ  homogenization  theory  such  as  composites.  This  is  because,  in  most 
cases,  the  ratio  of  scales  in  atomistic-continuum  problems  better  represents  the  asymptotic 
assumption  in  homogenization  than  in  other  types  of  multiscale  problems.  Furthermore,  the 
periodic  nature  of  nanopatterned  systems  is  more  apt  at  realistically  adhering  to  the  periodic 
assumption  than  unit  cell  models  of  composite  material  inclusions. 

In  this  report,  the  progress  and  initial  results  for  a  computational  framework  for  homoge¬ 
nization  of  the  atomistic  problem  is  presented.  Using  two  concurrent  domains,  one  for  the 
macroscale  continuum  domain  and  one  for  the  periodic  atomic  scale  domain,  self-consistent 
sets  of  equations  are  derived.  Atoms  in  arbitrary  configurations  and  structures  of  unlimited 
size  are  permitted.  Through  the  asymptotic  expansion  homogenization  technique,  a  set  of 
hierarchical  equations  are  derived  based  on  hyperelasticity.  At  the  local  level,  the  atomistic 
equations  are  used  under  the  assumption  of  the  harmonic  approximation  to  generate  the 
effective  properties  needed  to  solve  the  effective  global  level  equations.  The  Cauchy-Born 
rule  [13]  is  applied  to  the  atoms  to  enforce  the  gross  deformation  of  the  continuum  on  the 
atoms.  This  circumvents  the  need  to  apply  kinematic  constraints  by  making  use  of  the  weak 
averaging  properties  of  homogenization. 
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The  outline  of  this  report  is  as  follows.  In  section  2,  the  conventional  continuum  equations  are 
shown,  eventually  leading  to  a  variational  form  based  on  the  principle  of  virtual  work.  Then, 
in  section  3,  the  multiscale  equations  are  developed,  resulting  in  two  sets  of  equations  that 
govern  the  local  and  global  length  scales.  By  introducing  the  atomistic  potential  in  sections 
4  and  5,  the  details  of  the  atomistic  formulations  are  presented  and  cast  in  a  variational 
form  for  use  in  the  multiscale  homogenization  method.  In  section  6,  the  derivatives  of  the 
atomistic  energy  potential  needed  to  complete  the  derivation  of  the  method  are  provided  in  a 
general  form.  In  section  7,  one-  and  two-dimensional  (1-D  and  2-D)  demonstrative  examples 
are  shown.  Closing  remarks  are  discussed  in  section  8.  Additional  details  of  the  analytical 
derivatives  of  the  Tersoff-Brenner  Type  II  potential  were  presented  in  [14]. 


2.  Continuum  Formulations 


This  section  describes  the  kinematics,  stress  definitions,  and  linear  momentum  conservation 
laws  needed  to  develop  the  homogenization  method  from  atomistic  principles. 


2.1  Kinematics 


Consider  an  open  set  V  in  5ft3  that  deforms  to  the  configuration  v  in  5ft3.  Points  in  V  are 
denoted  X  =  (Xi,X2,Xs)  G  V  and  are  called  material  points,  whereas  points  in  v  are 
denoted  x  =  (aq,  x2, X3)  G  v  and  are  called  spatial  points.  The  deformation  is  a  one-to-one 
mapping  through  <f)  so  that  x  —  4>{X).  The  deformation  gradient  is  defined  by 


F  = 


d<p 

dx 


<9x  „ 

ax  =  VoX 


and 


d<f>i  _  dxi 
dX~j  ~  ~dXi 


(1) 


where  Vo  signifies  the  gradient  taken  with  respect  to  V.  The  determinant  of  F  is  termed 
the  Jacobian  and  is  defined  by  J  =  detF.  The  right  Cauchy-Green  strain  tensor  is  defined 
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by 


C  =  FtF,  (2) 

and  the  Green  strain  tensor  is  defined  by 

E  =  i(c-I).  (3) 


2.2  Stress  and  Equilibrium 


The  material  representation  for  the  conservation  of  linear  momentum  is  defined  by 

Vo-P  +  fo  =  0,  (4) 

where  P  is  the  first  Piola-Kirchoff  stress  tensor  and  fo  is  the  body  force  per  unit  of  undeformed 
volume.  In  rate  form,  it  is  given  by 


Vo  •  P  +  fo  =0. 

Using  the  principle  of  virtual  work,  equation  (5)  can  be  rewritten  as 

J  (v0  •  p)  6udV  +  J  f0  •  SndV  =  0,  V5u, 


(5) 


(6) 


where  5u  is  the  virtual  displacement.  Then,  using  the  definition  for  traction  with  respect  to 
the  undeformed  body,  equation  (6)  can  be  rewritten  as 


/  P  :  V0SudV  =  f  t0  •  6udA  +  f  f0 
Jv  Jav  Jv 


6udV. 


(7) 


We  invoke  the  notion  of  hyperelasticity  by  assuming  that  the  atomistic  potential,  W,  which 
is  a  function  of  the  atom  positions,  can  be  expressed  in  terms  of  strain.  This  assumes  that  the 
strain  energy  density  (or  the  free  energy  at  zero  temperature)  is  equivalent  to  the  atomistic 
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energy  potential.  Following  classical  continuum  mechanics,  one  can  then  define  the  first 
Piola-Kirchoff  stress  as 


P  = 


dW 

dF 


and 


p  =  SW 

v  SFij’ 


and  the  first  Lagrangian  elasticity  tensor  [15]  as 

d2W  _  dF 
“  dF 


C  = 


and  Cijki  — 


d2W 


dPij 


(8) 


(9) 


<9F<9F  dF  "  dFi:jdFki  dFkl 

A  relationship  is  needed  between  stress  and  strain.  From  equation  (9),  one  can  see  that  in 
hyperelastic  materials,  P  is  related  to  F  through 


P  =  CF  and  PtJ  =  CijkiFki, 


(10) 


where 

F  =  dii/dX  =  dv/dX,  (11) 

and  where  u  =  v  denotes  the  velocity. 


Substituting  equation  (10)  into  (7)  and  using  (11)  yields 

[  C  ::  (V0<5u®  V0v)dV=  f  t0-SudA+  [  f0  •  SudV,  V5u,  (12) 

Jv  JdV  JV 

and  the  equivalent  indicial  form, 

[  C‘i“^T7rrdV=  f  ioMdA+  [  foMidV.  (13) 

Jv  OAj  OAi  J dV  JV 

This  is  the  virtual  work  equation  associated  with  hyperelasticity.  The  two-scale  approach  is 
described  next.  It  is  devised  so  that  traditional  FE  continuum  equations  can  be  solved  in 
the  coarse  scale  and  atomistic  equations  can  be  solved  in  the  fine  scale. 


3.  Homogenization 


The  homogenization  framework  enables  the  weak  coupling  of  the  continuum  to  the  atoms. 
By  taking  the  limit  of  the  time-independent  asymptotic  expansion  parameter  e  — »  0,  we 
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exploit  the  weak  convergence  properties  of  the  scheme  in  order  to  decouple  the  length  scales. 
At  the  fine  scale,  the  domain  contains  only  atoms  with  periodic  conditions  prescribed  on 
the  boundary,  and  all  atom  displacements  are  measured  relative  to  a  fixed  point  in  the 
local  frame  of  reference.  From  classical  examples  of  continuum  mechanics  of  composite 
materials  [16],  this  enables  the  method  to  account  for  mutual  interactions  of  periodically 
spaced  heterogeneities  or,  in  this  case,  periodic  lattice  defects. 


The  homogenization  method  is  based  on  the  assumption  that  two  scales  exist  —  a  coarse 
scale  and  a  fine  scale.  Coordinates  in  the  coarse  material  scale  are  X  =  (XUX2,XZ),  and 
those  in  the  fine  material  scale  are  Y  =  Likewise,  the  spatial  coordinates  are 

the  lowercase  analogues.  The  two  scales  are  related  by  the  scale  parameter 

v  x 

y  =  7  <14> 

Therefore,  we  assume  that  the  ratio  of  scales  remains  the  same  before  and  after  deformation. 
The  aim  is  to  obtain  two  sets  of  coupled  equations.  The  asymptotic  series  assumption 
decomposes  the  displacements  as 


u(X)  =  U[°](X)  +  ut1)(X) 

=  ut0l(X)+euW(Y), 


(15) 

(16) 


where  represents  the  displacement  at  the  coarse  scale  and  uW  represents  the  perturbed 
displacements  due  to  inhomogeneity  at  the  fine  scale.  Square  brackets  denote  the  order  of  the 
term  in  the  asymptotic  series.  The  actual  physical  representation  of  the  total  displacement 
at  the  fine  scale  is  given  by  Takano  et  al.  [17]  as 


^u(X)  =  umicr°(Y) 

=  F(uM(X))Y  +  u[1](Y). 


(17) 


The  variable  X  in  equation  (17)  is  a  fixed  value  with  respect  to  Y.  That  is,  the  deformation 
gradient  of  a  point  in  the  coarse  scale  gets  mapped  onto  a  fine  scale  grid.  This  point  is 
typically  a  quadrature  point  in  an  FE  sense. 
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The  time  derivatives  are  analogous  to  equations  (16)  and  (17).  They  are  given  as 

u(X)  =  v(X) 

=  vl°l(X)  +  evU(Y),  (18) 

umicro(Y)  =  vmicr0(X) 

=  F(vl°l(X))Y  +  vW(Y).  (19) 


(20) 


(21) 


Substituting  equations  (16)  and  (18)  into  (13)  yields 

[  C  ::  [Vx  (<5u[0](X)+  £<5uW(Y))  (v^X)  +  evW(Y))]  dV 

Jv 

=  f  (<5u[0]  (X)  +  e<5u[1]  (Y))  •  t0cL4 
JdV 

+  [  (W°]  (X)  +  eSu[1]  (Y))  •  f0 dV,  VW0] ,  <5u[1] . 

Jv 

Note  that  by  use  of  the  chain  rule  and  equation  (14), 

<9Y 

Vx0(X,Y)  =  Vx0+^Vy0 

=  H - Vy</>. 

£ 

Therefore, 

V*  (u[0](X)  +  eu[1](Y))  =  V^u[°](X)  +  VyuW(Y). 

Using  equation  (22)  in  (20)  and  taking  the  average  over  Y  gives 

[  tL  [  C::  [(Vx<5u'°](X)  +  VySu^(Y))  0  (Vxv^(X)  +  VyvW(Y))]  dYdV 
Jv  1*1  Jy 

=  [  (W°]  (X)  +  e<5u[1]  (Y))  •  t odA 
Jav 

+  f  (<5u[0]  (X)  +  eSu[1]  ( Y) )  •  f0dV,  V<5u[0] ,  Su[1] . 

Jv 

Then,  in  the  limit  as  e  — >  0,  equation  (23)  is  satisfied  only  if  the  following  two  equations  are 
satisfied, 

[  [  C  ::  [Vx«5uM(X)  ®  (VxvM(X)  +  VyvW(Y))]  dYdV 
I*  I  Jv  Jy 

=  [  <5u[0](X)  -i0dA  +  [  <5u[0](X)-f0dU,  V  <5ul°], 

JdV  Jv 


(22) 


(23) 


(24) 
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W\LLC  ::  [Vy<5utll(Y)  ®  (V^vt0](x)  +  VyvW(Y))]  dYdV  =  0,  V  <5uW.  (25) 

By  recourse  to  the  FE  method,  the  solution  of  equation  (24)  is  straightforward,  assuming  C 
and  vM  are  known.  It  is  then  evident  that  due  to  the  dependence  of  equation  (25)  on  v^, 
equations  (24)  and  (25)  are  coupled  and  must  be  solved  concurrently.  For  general  problems, 
an  iterative  numerical  solution  scheme  can  be  employed  to  handle  the  non-linear  system  of 
equations  together  with  a  linearly  ramped  load  to  ensure  solution  convergence. 

In  the  next  section,  a  method  is  shown  for  solving  equation  (25)  for  v[1l  Then  in  the  following 
section,  the  formulation  that  enables  the  atomistic  information  to  be  fed  into  equation  (24) 
is  derived.  These  two  sections  constitute  the  iteration  steps  that  must  be  performed  for  a 
general  application. 


4.  Atomistic  Equation 


Distinct  and  distinguishable  atoms  are  assumed  to  reside  in  the  local  level  cell.  By  the 
Cauchy-Born  rule  [13],  at  a  point  X,  F(u^)  is  assumed  to  give  the  energy  minimizing 
configuration  of  the  atoms.  For  simplicity,  we  assume  that  the  atoms  are  arranged  in  a 
lattice.  Then,  the  positions  of  the  atoms  Y  are  given  from  the  lattice  coordinates  m  by 

Y(m)  =  me,  :  m  G  £,  £  =  Z3,  Z  <  N,  (26) 

where  e*  are  the  primitive  translation  vectors  and  N  is  the  integer  multiple  of  atoms  con¬ 
tained  in  the  unit  cell.  To  avoid  confusion  in  notation,  atom  labels  are  noted  in  parentheses 
henceforth  and  are  not  subject  to  the  conventional  summation  rules  associated  with  indicial 
notation.  The  displacement  of  the  atoms  are 

q(m)  :  m  6  £.  (27) 

Note  that  there  is  no  restriction  to  perfect  lattices.  In  fact,  by  using  computers,  arbitrary  arrangements 
of  atoms  can  be  considered  as  long  as  the  assumption  of  the  Cauchy-Born  rule  still  applies. 
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Upon  deformation,  the  new  positions  of  the  atoms  are  given  by 

y(m)  Y(m)  “b  Q(m)- 

The  deformation  gradient  is  defined  by 


(28) 


(29) 


The  vector  separating  two  atoms,  i  and  j,  in  the  reference  configuration  is  given  by 

R(y)  =  Y(j)  —  Y(j), 


(30) 


where  Yq-)  denotes  the  position  of  atom  j\  and  Y(q,  the  position  of  atom  i.  The  vector 
separating  two  atoms  in  the  deformed  configuration  is  given  by 

r(ij)  =  y u)  -  y «■  (31) 


Then,  the  Cauchy-Born  rule  can  be  stated  in  a  more  precise  manner  by 

r(u)  =  — 

=  FR(ij). 

In  defect  regions  and  through  the  homogenization  theory  via  equation  (17),  the  rule  becomes 


(32) 


r(ij)  FR(ij)  T  , 

where  f^j)  =  -  u^j  is  the  additional  term  to  account  for  high  energy  regions. 


(33) 


For  the  energy  associated  with  the  deformation  of  the  atoms,  we  use  a  modified  form  of  the 
so-called  Potential  II  parameterization  of  the  Tersoff-Brenner  potential  [18,19].  It  takes  the 
form 

W=-b[Bt(Y  +  q)-Bl(Y)],  (34) 
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where  W  is  the  energy  density  of  the  frozen  system,  N  is  the  number  of  atoms,  and  Eb  is 
the  binding  energy  given  for  a  pure  carbon  system  by 


=  EEWy-iM'i#))]. 

i  j(>i) 

B  —  2  (^(u)  **’  B(ji))  , 

v(r)  =  /w)(r)Pw  -%/2S/?(r-.R('>) 

K{'  (S- 1) 

|/4(r)  =  fm(r)D{e)S 
Ml  (S-l) 


(35) 

(36) 

(37) 

(38) 


/(«)(’■) 

EUJ) 

cm 


1,  r  < 

<  i{l  +  cos[fJ^]},  Rm<r<  RW 

k  o,  r  > 

-1  -<5 

1  +  E  G(0(ijk))f(ik)  (r(ik)) 

a  /i  +  4  *  } 

+  d2  d2  +  (l+cos0)2/’ 


(39) 

(40) 

(41) 


with  the  constants  given  in  Table  1 .  The  modification  in  this  work  comes  by  way  of  omitting 
the  extra  bond-order  term  in  [19],  which  is  primarily  used  for  problems  involving  changes 
in  coordination  numbers.  We  therefore  presently  restrict  our  consideration  to  classes  of 


deformation  involving  no  change  in  coordination. 


Given  that  the  energy  can  be  written  as  a  function  of  the  atom  displacements,  equation  (25) 
can  be  expressed  in  a  form  conducive  to  atom  representations.  We  equate  vM  to  the  rate  of 
atom  displacement  and  attempt  to  solve  the  equivalent  form 


d_n  c>41]  _  acmd 40] 
dY3  l]kl  dYt  dY3  dXx 


(42) 


under  periodic  boundary  conditions.  The  solution  to  equation  (42)  is  found  as  the  zero  of 
&R,  in  the  equation 


dn  =  TCvW  -  V  •  V0v[°J(x), 


(43) 
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Table  1.  Parameters  for  Tersoff-Brenner  potential. 


1.39  A 

£)(e) 

6.0  eV 

S 

1.22 

(3 

2.1  A 

S 

0.5 

RW 

1.7  A 

rw 

2.0  A 

cto 

0.00020813 

cl 

3302 

dl 

3.52 

where  K  is  the  N  x  N  Hessian  and  is  given  by 

d2W 

K  =  7—-,  (44) 

dqdq 

where  q  is  the  vector  of  atom  displacements  of  size  3 N  (in  three  dimensions)  and  X>  is  a  third- 
order  unsymmetric  tensor  that  is  obtained  from  the  first  derivative  of  the  Euler-Lagrange 
equation  with  respect  to  the  local  deformation  gradient  given  by 

_  d2W 

V~  dqdF  ^ 

The  size  of  T>  depends  on  the  dimensionality  of  the  problem.  In  three  dimensions,  it  can  be 
expressed  as  an  N  x  9  matrix,  where  9  corresponds  to  the  number  of  independent  components 
of  F. 


5.  Multiscale  Equation 


Once  equation  (43)  has  been  solved  for  vM,  the  remaining  task  is  to  formulate  a  tractable 
global  scale  boundary  value  problem.  The  key  distinction  between  this  investigation  and 
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conventional  continuum  formulations,  such  as  hyperelasticity,  is  the  conspicuous  incorpora¬ 
tion  of  vM,  a  fine-scale/atomistic  quantity,  in  the  global  scale  equations,  and  the  definition 
of  the  material  property  tensor  completely  in  terms  of  atomistic  variables. 


We  return  to  equation  (24),  recognizing  that  vW  is  now  known.  Incorporating  the  definition 
for  the  first  Lagrangian  elasticity  tensor  from  equation  (9)  yields 

WlUyWd f  I7*'5""’'®  ®  (V*vM(X)  +  Vvvl'l(Y))]  dYdV 


[  <5u[°](X)  •  1 0dA  +  f  <5u[0](X)  •  f0dV,  V  6u[0l 
Jav  Jv 


Then,  using  the  definition  of  F  in  equation  (1)  and  assuming  a  first-order  Taylor  series 
representation  for  the  time  derivative  gives, 

Iff  d2W 

m  L  L  otw  ::  ®  VxvW(X))  dYdV 


dFdF 


=  £W°1(X)-M4  +  J  W°](X)  -f0dV  (47) 

Iff  d2W 

~W\JvJy8m:'  (V^'5u""(X)  ®  vll|<Y)) dYdV-  V  ■S“W. 

where  there  is  now  a  double  contraction  on  and  a  single  contraction  on  vT  jn  the 

last  expression  of  equation  (47).  The  solution  to  equation  (47)  yields  vt°l  It  is  noteworthy 
that  the  last  term  is  zero  when  the  energy  distribution  over  Y  is  constant,  i.e.,  when  the 
atom  arrangement  forms  a  perfect  lattice.  This  reduces  the  problem  to  a  classical  harmonic 
approximation  where  the  first  Lagrangian  elasticity  tensor  is  assumed  to  model  the  material 
behavior.  In  equation  (47),  the  last  term  serves  as  a  corrective  force  in  regions  of  highly 
energetic  atoms  (i.e.,  nonlocal  regions)  to  account  for  defects  and  lattice  inhomogeneities. 


•fo  dV 


:  ■  (V*<5ut°](X)  (2)  vW(Y))  dYdV,  V  <5u 


6.  The  Euler- Lagrange  Equations  and  the  Hessian 

In  this  section,  the  analytic  forms  of  the  Euler-Lagrange  equations  and  Hessian  are  derived 
for  a  general  potential.  The  nontrivial  algebra  typically  needed  to  obtain  equations  (44)  and 
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(45)  for  the  specific  case  of  the  Tersoff-Brenner  potential  are  shown  in  greater  detail  in  [14], 
and  only  general  forms  are  derived  here.  The  Euler-Lagrange  equation  is  the  first  derivative 
of  the  Lagrangian  with  respect  to  the  degrees  of  freedom.  In  this  problem,  the  Lagrangian 
is  the  negative  of  the  atomistic  energy  density,  which  we  presently  assume  is  equivalent  to 


the  free  energy  at  zero  temperature.  The  Euler-Lagrange  equation  is  therefore  given  by 

dW 


£  = 


dq  (m) 

1  dEb 


(48) 


N  dq (m)  ’ 


and  using  the  chain  rule  for  derivatives,  it  is 

1  dEb 


£  = 


N  dq (TO) 
dEb 


i( 


dr 


(b) 


+ 


dEi 


L  df(it)  +  aE*  .  8r(jt)^ 


(49) 


(jj)  ^Q(m)  &£(ik)  ^Q(m)  <9Q(m)  / 


Here,  we  have  implicitly  assumed  that  there  are  three  independent  atomic  position  vectors. 


One  can  show  quite  easily  that  there  are  in  fact  only  two  by  using  the  relationship  ryfc)  = 


r(ik)  r(b)  ■ 


The  Hessian  is  obtained  by  taking  an  additional  derivative  of  the  Euler-Lagrange  equations. 

Specifically,  we  again  make  use  of  the  chain  rule  to  obtain 
d2W 


K  = 


dq(n)dq(m) 

1  d^Eb 
N  dq(n)<9q(m) 
d2Eh 


N  [dr(ij)dr{ij) 
d2Eb 


dr(ij)  q  ^r(b) 


+ 


+ 


+ 


dr  (itydr  (ij) 
d2Eb 


dq(n) 

( 

V&Kn) 

.  /  dr(3k) 


® 


<g> 


dq(m) 

dq(m)J 


d2Eh 


dr 


(b) 


dr^jk)dr{ij) 

d2Eb 


\  ^0.(91)  ^Q(m) 


+ 


d2Eh 


+ 


dr(jk)dr(jk)  '  \d q(n)  dq(m)  J 


+ 


dT(ij)d*(ik) 

dr{ik)dr{ik)  ' 
d2Eb 


( 


dr(ij)  ^  dr(ik)  \ 

dq(m)  J 


+ 


dr{ij)dr(jk) 

d2Eb 


dq{n) 
(  drm 
V^q(n) 
Edr^i 
\dq(n) 
dr(ik) 


® 


® 


’( 


dq(m)J 

dq(m) 
d T(jk ) 


dr{ik)dr(jk)  '  \^q(n)  9% 


drW  Q  dr{lk)  \ 
dr{jk)dr{ik)  '  \  dq{n)  dq(m)  J  J 


(50) 
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Second  derivatives  of  the  interatom  vectors  are  zero,  i.e., 
d2r{ij)  d2 


dq(m)dq(n) 

=  0  V(m,  n). 


(Y(»)  +  qo)-Y(o  +  q(i)) 


(51) 

(52) 


Next,  the  appropriate  right-hand  side  expressions  are  derived  for  equations  (43) 
This  involves  the  use  of  the  chain  rule  again  to  obtain 


and  (45). 


<92W 

dq(ra)dF 


N 


■( 

dr(lk)dr{jk)  '  ( 


dr{ij)dr{lj) 
,  &Eb 
dr(ik)d*(ij) 
d2Eb 

dT(ij)dT{jk) 

d2Eb 


drfa)  0  &r(v)  ^  d2Eb 


OF 

dF 

dF 

dF 


+ 


dq  (rn)J  dr(ij)dr{ik) 

0  )  |  ^ 
^q(m)/  dr^dr^ 

5rM\  _  d2£6 


and  by  definition, 


drfa) 

dF 


'j  and 


+ 


dq  (TO)  / 


d2Eb 


d(l  (m)J  dr{jk)dr{ik) 

,  d2£„ 


dr^ndr 


(jk)ar(jk) 


di(il 


m 


dF 


=  R, 


(»*)  • 


(53) 


(54) 


Finally,  we  use  a  similar  approach  to  define  the  first  Lagrangian  elasticity  tensor.  This  is 

the  traditional  way  of  estimating  the  elastic  properties  of  a  solid.  Using  the  chain  rule  once 
again  gives 

d2W 
dFdF 


d2Eh 


dii 


AT 

Ldr(y)dr{ii)  ‘ 

lap 

•  ® 

ap  ) 

4. 

d2£6 

/"  8r<(fl 

® 

T 

di(ik')di^{j'j 

V  dF 

dF  ) 

_L . 

d2Eb 

® 

di{ij)\ 

dr(ij)di(jk) 

V  dF 

dF  J 

J - 

d2Eb 

(  drUk) 

0 

dr(ijt)  \ 

dr(ik)dT(jk)  ‘ 

V  dF 

dF  J 

d2  Eh 


+ 


+ 


di  (ij)dx  (^ik^ 

d2Eb 

di{ik)dim 

d2Eb 

^v(jk)dl{ij) 

d2Eb 

dr^^di^ 
d2Eb 


_ .  ( dr(*U  ~  d*(v)  \ 

ik)  '  V  dF  ®  dF  ) 


dF 

( g£w  drm\ 
\  dF  ®  dF  ) 


dF 

di(ik) 


( 


dF 

d£M 

dF 

d*(ik) 

dF 


® 


® 


dF 

dF 


9rW5rW  '  (  dF  ™  dF 


dF 


) 


(55) 
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The  first  Lagrangian  elasticity  tensor  is  used  in  equation  (47),  whose  solution  gives  vt°l  In 
a  perfect  lattice,  equation  (55)  provides  the  only  atomistic  material  information  needed  to 
solve  the  macroscopic  continuum  problem.  The  next  section  illustrates  this  by  showing  that 
the  perturbation  is  zero  for  a  uniform  crystal. 


7.  Example  Problems 

7.1  Example  I:  Perfect  1-D  Atomic  Lattice 

To  illustrate  the  calculation,  a  1-D  analytical  example  is  presented.  The  Tersoff-Brenner 
potential  is  used  to  represent  the  energetics  of  a  1-D  single-species  chain  of  carbon  atoms. 
The  objective  here  is  to  solve  equations  (42)  and  (43)  for  and  demonstrate  a  simple  case 
of  a  perfect  lattice  using  this  method. 

One  atom  comprises  the  periodic  unit  cell,  but  to  account  for  the  effects  of  triples,  two 
“fictitious”  atoms  are  assumed  to  extend  beyond  the  boundaries  of  the  cell  on  each  side  as 
illustrated  in  Figure  3.  Periodic  conditions  apply  at  the  cell  boundaries.  The  equilibrium 
lattice  constant  for  the  chain  is  rQ  =  1.86868A. 


Fictitious  atoms 
/  \ 

Fictitious  atoms 
/  \ 

O 

O 

i  •  i 

O 

O 

5 

3 

:  1  : 

\  /' 
Cell  boundaries 

2 

4 

Figure  3.  Unit  cell  of  the  1-D  carbon  chain.  The  atoms  are  labeled  by  iden¬ 
tifying  numbers. 
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The  following  two  conditions  stem  from  the  1-D  assumption, 


6  =  7T, 

RW  <  r  <  RW. 


(56) 


This  simplifies  the  expressions  in  [14].  The  resulting  Hessian  for  three  arbitrary  colinear 
atoms  ( i ,  j,  k )  is  obtained  as 


[K{ijk)]  = 


y{ijk)  y~{ijk)  ijk ) 

Ml  M2  M3 

v-(ijk)  v~(ijk)  y{ijk ) 

Ml  M2  M3 

v-(ijk)  y{ijk)  y{ijk) 

/Vqi  A^32  /V/qq 


'33 


(57) 


where  JCmfi  =  JCnffl  and  the  terms  are  defined  by 

>  =  VR-  BV';  +  +j«+ 


dpt  f" 

+  2  ^ 

tAijk)  _  y"  Dy"  D1+J  f' 

^12  “  “  ~YVAB(ij)  ! [iky 

(ln& 


K“*,_  2  V'|Jt!)A«)"5(f  + W(S‘  («4ii)  ~  ■yv'-4B(.j|I/wi> 


4f  ’  =  vR  -  gy;, 

/r(b'fc)  _  ao^  t/  d1+j  f' 
^23  -  ~2~VAB(ij)  j(ik)i 


4f 1  =  *(«  +  dk,b$  (a./w)2  +  fvAB$f;ty 


(58) 

(59) 

(60) 
(61) 
(62) 
(63) 


Upon  assembly  of  the  two  unique  pairs  (31,12)  and  their  associated  triples  (123,132)  in 
Figure  3,  the  final  assembled  Hessian  of  the  global  system  is  given  by  the  matrix, 


K,\\  /C12  /C13 

=  K-21  /C22  /C23 

^31  ^32  ^33 

which  is  assembled  through  the  operation, 

W  =  U  U  I^*’]  =  Km«, 

m)  (n) 


(64) 


(65) 
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where  U  is  the  addition  operator  over  all  unique  pair  and  triple  combinations  of  (i,  j,  k)  and 
(m)  and  (n)  are  displacement  degrees  of  freedom  for  each  atom.  In  equation  (65),  [K]  is 
symmetric  once  again,  and  its  components  are  obtained  in  detail  for  the  problem  shown  in 
Figure  3  as  follows: 


Kn  =  Kg2”  +  +  Kgil\ 

(66) 

K12  =  Kg23’  +  Kg“>, 

(67) 

y  i^(123)  .  k(315) 

A-13  —  A^13  +  1^12  , 

(68) 

K*  =  Kg23’  +  Kg41’. 

(69) 

*^(123) 

^23  —  A^23  j 

(70) 

K33  =  4f)+Ki3115). 

(71) 

This  constitutes  the  stiffness  matrix  K  in  equation  (43). 


The  next  step  is  to  calculate  the  right-hand  side  of  equation  (42),  which  is  equivalent  to  cal¬ 
culating  T>  and  multiplying  by  the  global  rate  of  the  deformation  gradient.  Using  equations 
(45)  and  (53),  the  right-hand  side  for  three  colinear  atoms  ( i,j,k )  is  obtained  as 


{V{ijk)}  =  i 


VW) 


\ 


n 


(ijk) 


(72) 


where  the  components  are  defined  by 

©P  =%>  (Vr  -bv;)  +  (%,,  -  fi(w) 

+  R{ik)  (2  (*5  +  1)VaB^s  (tto/(ife))  +  /(**))  . 

=  -  Riij)  {K  -  Bvj)  -  R(ik)  (j?-Vab$  /,'*,)  , 

=  -  Rrn )  (x^B<w  ^)) 

-  Rm  (|(«  + 1  (“»/w)2  +  °4V* ■ 


(73) 

(74) 

(75) 
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As  earlier,  the  assembly  operation 


(hj,k) 

{D}=  LI  {!>«*)}  =D„, 

(m) 

yields  the  right-hand  side  of  the  global  system  given  by, 

M 

{©}  =  <  v2  L 

where  the  components  are 


(76) 


(77) 


©,=I><m)+X>fs>  +  P«241',  (78) 

V2  =  V"n>  +  pj21’>,  (79) 

©a  =  2><m>  +  ©<”>.  (80) 

Under  the  assumption  of  a  1-D  perfect  lattice,  we  have  and,  consequently, 

Vi  =  0.  Then,  we  can  satisfy  the  periodicity  condition  and  the  rigid  body  constraint  by 
setting  =  0.  The  solution  is  therefore 


VW  -  „W  - 

v(i)  ~  V(2)  ~  V(3) 


0. 


(81) 


In  light  of  equation  (81),  the  last  term  in  equation  (47)  is  zero,  and  the  material  properties  are 
obtained  from  the  atomistic  energy  density  solely  through  equation  (55).  This  result  shows 
that  in  a  defect-free  lattice,  the  homogenization  method  coincides  with  the  conventional 
atomistic  hyperelasticity  problem.  The  next  section  shows  an  example  in  which  a  defect 
causes  an  inhomogeneous  energy  distribution,  leading  to  a  situation  where  homogenization 
is  needed  to  average  out  the  energy. 


7.2  Example  II:  1-D  Atomic  Lattice  With  Defect 


Consider  the  problem  shown  in  Figure  4,  where  the  center  atom  is  displaced  by  a  distance 
L  from  its  original  energy  minimizing  configuration.  This  displacement  of  the  center  atom 
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Fictitious  atoms 

Fictitious  atoms 

/ 

\ 

L  . 

/ 

\ 
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O 

O#  i 

O 

GO 
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3 

1  : 

2 

4 

\ 

/' 

Cell  boundaries 

Figure  4.  Unit  cell  of  one-dimensional  carbon  chain  with  periodic  defect, 
constitutes  the  defect. 


With  this  change,  the  key  stiffness  matrix  term  in  equation  (66)  is  now 


and  likewise,  equation  (78)  is  now 

#1  =#(12)  (y^12)  -  #(  12)K4(12))  +  (#(13)  -  #(12))  (^YVAa2)B}u)f(  13)) 

+  #(13)  Q(<5  +  l)^A(12)#(1it)<  (G°/(  13))  +  ~yVa(12)B(^  f{  13)^ 

-  #(13)  (Vr(13)  -  #(13)  V^(13))  -  (#(12)  -  #(13))  ^~|"^4(13)#(13)6/(12)^ 

-  #(12)  Q(<5  +  1)Va(13)#(1+3)  (a°/(  12))  +  ~YVaH 3)#(13*  /(12))  > 


(82) 


(83) 


where  #(i2)  =  rG— L  and  #(13)  =  r0+L.  Then,  solving  equation  (43)  under  periodic  boundary 
conditions  gives 


[X]  =  gtg 
1,1  /Cx  5a; 


(84) 


The  solution  as  a  function  of  Ljr0  is  shown  in  Figure  5.  As  expected,  the  solution 

has  symmetry  about  the  origin  and  grows  asymptotically  larger  as  the  size  of  the  defect  (L) 
grows  closer  to  the  cut-off  radii.  We  intentionally  avoid  larger  defects  due  to  the  nonconvex 
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structure  of  the  energy  well  associated  with  the  Tersoff-Brenner  potential.  This  generally 
leads  to  unphysical  discontinuities  in  the  perturbation  velocity  (v^)  due  to  discontinuous 
second  derivatives  of  the  atomistic  energy  with  respect  to  the  defect  size.  This  is  attributable 
to  the  construction  of  the  empirical  potential  in  equations  (35) — (41 ) ,  which  is  suited,  by 
design,  for  systems  where  nearest  neighbor  atoms,  even  in  defect  regions,  are  within  the 
cut-off  radius 


Figure  5.  Distribution  of  vW/VqvW  solution  as  a  function  of  the  defect  size 
L. 


It  is  also  noteworthy  that  arbitrary  defect  densities  can  be  treated  by  appropriate  modifica¬ 
tion  of  the  unit  cell.  In  most  cases,  one  can  tailor  the  desired  density  by  increasing  the  size 
of  the  unit  cell  and  performing  the  summations  and  the  assembly  of  the  atomistic  discrete 
equations  over  more  atoms.  Figure  6  illustrates  this  idea  for  the  1-D  carbon  chain. 

Numerical  experiments  show  that  as  the  size  of  the  unit  cell  increases,  the  perturbative 
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Figure  6.  Larger  chain  of  atoms  in  perfect  arrangement  around  the  defect 
region  decreases  the  defect  density. 

displacement  has  a  sharp  discontinuity  at  the  defect.  Figure  7  shows  this  nonlocal  behavior 
as  the  number  of  atoms  increases.  The  problem  is  of  a  single  defect  in  chains  of  increasing 
size.  The  defect  magnitude  is  held  fixed  at  L/ra  =  0.01.  The  nonlocal  discontinuity  of  the 
perturbative  velocity  qualitatively  agrees  with  traditional  displacement  jumps  that  occur  at 
dislocation  cores.  The  discontinuity  indicates  that  the  material  property  at  the  defect  ( Jf^) 
is  modified  by  the  last  term  in  equation  (47),  an  amount  proportional  to  i;W  that  serves  as 
a  correcting  force  for  the  nonlocal  effect. 


Figure  7.  Distribution  of  uW/Vou^  along  unit  cell  length  for  varying  num¬ 
bers  of  atoms  ( L[r0  =  (J.U1). 
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Although  the  primary  details  of  the  method  have  been  demonstrated  in  these  two  examples, 
the  method  can  be  extended  to  consider  the  multiscale  problem  show'n  in  equation  (47) 
for  more  general  cases  involving  self-consistent  solutions  with  equation  (43)  as  in  the  next 
section. 


7.3  Example  III:  2-D  Graphene  With  Defect 

The  method  can  be  generalized  to  multiaxial  problems.  In  this  example,  we  consider 
graphene  with  three  types  of  point  defects:  interstitial,  equilibrium  point  vacancy  and  sad¬ 
dle  point  vacancy.  The  atom  positions  around  the  point  defects  were  first  computed  using 
quenched  molecular  dynamics.  These  are  illustrated  in  Figure  8  together  with  the  original 
defect-free  configuration.  Similar  defect  structures  have  been  encountered  both  numerically 
and  experimentally  [20-22]. 

The  effective  elastic  constants  for  the  defect-free  case  were  in  the  form  of  the  first  Lagrangian 
elasticity  tensor  (55),  which  compare  reasonably  with  experimental  results.  Using  the  second 
derivative  of  the  Tersoff-Brenner  potential,  the  bulk  values  (in  units  eV/atom)  for  graphene 
were  computed, 


Cmi  =  66.51 

C2112  =  21.63, 

(85) 

C1122  =  20.06 

C2121  =  24.83, 

(86) 

Cun  =  24.83 

£-221.1  =  20.06, 

(87) 

C1221  =  21.63 

£2222  =  66.51, 

(88) 

and  terms  not  listed  are  zero.  The  equilibrium  energy  is  -7.37563  eV/atom  and  nearest 
neighbor  bond  length  1.45  A.  For  an  assumed  layer  thickness  of  3.4  A,  which  is  the  standard 
layer  separation  thickness  for  graphite,  the  effective  Young’s  modulus  from  the  bulk  is  Y  = 
1.261  TPa  with  an  effective  Poisson’s  ratio  of  0.302.  These  values  agree  well  with  measured 
values  for  graphite  and  carbon  nanotubes  [23]. 
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Figure  8.  Atom  configurations  in  periodic  cell.  Lines  are  used  to  denote 
regions  near  the  defect. 

The  model  problem  solved  at  the  macroscopic  scale  is  depicted  in  Figure  9.  The  uniform 
grid  is  composed  of  25  4-noded  quadrilateral  elements.  The  right  edge  is  pulled  uniformly 
and  the  left  edge  is  held  fixed.  All  units  of  measure  are  carried  through  in  terms  of  eV  and 
A  so  that  no  assumption  of  a  layer  thickness  is  required.  That  is,  although  the  problem  is 
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Figure  9.  FE  model  of  macroscopic  problem. 


two-dimensional,  there  is  no  need  for  a  plane  assumption. 

The  solution  procedure  involves  application  of  incremental  loads  with  iterative  loops  over 
equation  (47)  and  sub-loops  over  (43).  The  sub-loops  seek  converged  values  for  vM  subject 
to  the  modification  from  equation  (33).  In  this  problem,  an  overly  cautious  load  increment 
of  1.0  x  10~4  strain  was  used,  which  ensured  convergence  tolerances  for  both  and  vM  far 
beyond  10-8A. 

The  effective  strain  energy  density  and  material  properties  are  shown  vs.  strain  in  Figures  10 
and  11.  Of  key  interest  are  the  sudden  jumps  that  occur  for  the  effective  material  properties 
of  saddle  point  and  interstitial  defect  scenarios  at  approximate  strain  values  of  0.05  and  0.15, 
respectively.  This  is  attributable  to  the  relative  instability  of  those  types  of  point  defects. 
Such  observations  were  made  previously  for  the  saddle  point  vacancy  in  [20].  The  effective 
properties  and  energies  at  larger  strains  after  the  initial  instability  are  meaningless  and  are 
therefore  not  shown.  The  instability  occurs  in  this  problem  because  the  separation  distance 
between  neighboring  atoms  near  the  point  defect  exceeds  the  cut-off  radius  of  the  potential 
because  of  the  deformation,  thereby  distorting  the  strain  energy  and  effective  properties. 


We  expect  that  the  actual  strain  is  much  smaller  before  the  instability  occurs  because  of  the 
present  zero-temperature  assumption. 


(a)  Cun 


(b)  C2222 


Figure  10.  Material  property  change  with  applied  uniform  strain. 


It  is  noteworthy  that  the  approximate  atom  density  in  the  longitudinal  direction  decreases 
(cell  elongates),  while  it  increases  in  the  transverse  direction.  The  subsequent  material 
nonlinear  effect  causes  the  properties  to  decrease  in  Cun  while  increase  in  C2222-  We  also 
note  that  the  interstitial  point  defect  possesses  the  highest  values  for  material  properties 
and  lowest  strain  energy,  whereas  the  trends  for  the  saddle  point  vacancy  are  precisely  the 
reverse.  The  defect-free  structure  exhibits  stronger  material  properties  and  higher  strain 
energy  than  the  equilibrium  vacancy  structure. 

Comparing  convergence  rates  between  methodologies  with  homogenization  (i.e.,  vW  ^  0) 
and  without  homogenization  (i.e.,  vM  =  0)  shows  the  former  with  a  distinct  advantage.* 

*The  methodology  without  homogenization  implies  an  approach  whose  material  property  tensor  is  com¬ 
puted  directly  from  the  second  derivative  of  the  energy  potential  and  whose  formulation  does  not  involve  a 
perturbative  term. 
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Figure  11.  Strain  energy  with  applied  uniform  strain. 

The  equilibrium  vacancy  configuration  is  used  here.  The  convergence  results  are  illustrated 
in  Figure  12.  The  calculation  is  of  the  transverse  displacement  of  the  upper  right  corner  node 
of  the  mesh  in  Figure  9  as  it  evolves  with  the  total  number  of  computation  cycles,  Nsteps. 
The  computation  cycle  is  computed  by  adding  the  total  number  of  nonlinear  iterations  steps 
(major  loops  plus  sub-loops)  with  the  total  number  of  load  increments,  which  is  a  crude 
way  to  estimate  the  convergence  behavior.  The  number  of  load  increments  is  selected  to 
ensure  that  the  final  transverse  displacement  in  the  two  methods  is  less  than  1%  different  in 
magnitude.  Figure  12  shows  that  the  total  required  number  of  calculation  steps  is  smaller 
in  the  homogenization  result  by  a  factor  of  four. 

This  convergence  behavior  is  attributable  to  the  addition  of  the  vW  term  in  equation  (33), 
which  is  closer  to  the  energy  minimizing  configuration  than  equation  (32)  alone.  Figure 
13  illustrates  the  atom  displacement  due  to  equation  (32)  relative  to  a  local  coordinate 
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Figure  12.  Convergence  with  and  without  homogenization. 

system  centered  on  the  periodic  cell,  and  Figure  14  illustrates  the  subsequent  correction 
from  equation  (33). 

The  uniform  strain  applied  in  this  example  results  in  a  uniform  state  of  strain  throughout 
the  entire  mesh  in  Figure  9.  It  should  be  noted  that  more  complex  loading  scenarios  can 
be  treated  merely  by  changing  the  boundary  conditions  of  the  problem.  These  simple  re¬ 
sults,  however,  are  indicative  of  the  generality  of  the  method  for  two-  and  three-dimensional 
problems. 


8.  Conclusions 


Linking  atomic  scale  physics  with  continuum  scale  phenomena  is  of  keen  interest  in  the 
mechanical  study  of  solids  and  nanostructures.  The  effects  that  dominate  the  mechanical 
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Figure  14.  Atom  displacements  due  to  homogenization  method. 


behavior  at  the  continuum  scale  typically  initiate  and  evolve  from  the  atomic  scale.  Moreover, 
periodic  structures  can  emanate  from  nanopatterning  and  epitaxy  through  stresses  induced 
from  an  underlying  substrate.  Despite  numerous  promising  methods  in  the  literature  that  are 
capable  of  linking  scales  up  to  the  micron  level,  periodic  structures  with  global  dimensions 
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at  and  beyond  the  millimeter  range  -  needed  for  mass  producing  nanoscale  devices  -  have 
only  begun  to  be  studied.  To  this  end,  we  have  attempted  to  address  this  issue  by  exploiting 
the  features  of  homogenization  theory  to  devise  a  scheme  that  passes  atomistic  information 
to  very  large  continuum  scales. 

We  have  applied  the  Cauchy-Born  rule  to  the  atom  scale  by  assuming  that  the  configuration 
of  atoms  used  to  solve  for  the  perturbation  displacement  at  each  load  increment  is  indeed 
the  minimizing  configuration  of  the  atomistic  energy.  We  have  not  considered  the  method 
in  conjunction  with  a  lattice  statics  routine,  i.e.,  various  strategies  of  minimizing  the  atom¬ 
istic  energy  by  quenched  molecular  dynamics  or  solving  Newton’s  equations  to  minimize  the 
interatom  forces.  These  assumptions  (i.e.,  zero  temperature  and  Cauchy-Born)  also  pre¬ 
clude  the  formulation  from  modeling  thermally  activated  phenomena  such  as  crack  growth, 
propagation,  or  damage  evolution.  However,  it  can  be  used  to  estimate  mechanical  effects 
across  coupled  length  scales  and,  if  needed,  serve  as  the  underlying  framework  for  modified 
algorithms  that  can  account  for  such  problems. 

We  demonstrated  the  method  for  one-  and  two-dimensional  problems  containing  point  de¬ 
fects.  Mechanical  data  were  reported  from  numerical  experiments.  The  paucity  of  exper¬ 
imental  data  for  nanomechanics  makes  validation  difficult.  However,  the  material  proper¬ 
ties  stemming  from  the  reference  configuration  compares  very  well  with  available  published 
measurements  and  the  observed  trends  from  mechanical  deformation  agree  with  generally 
accepted  intuition. 

For  this  work,  the  specific  case  of  the  Tersoff-Brenner  Type  II  potential  was  considered. 
But  the  principles  and  the  general  equations  can  be  extended  to  any  potential,  provided  the 
appropriate  derivatives  can  be  obtained  as  in  [14].  Typically  for  classical  systems,  onerous 
tensor  algebra  and  calculus  are  required  or  more  computationally  efficient  procedures  can  be 
implemented  to  obtain  derivatives  numerically  [24]. 
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The  aim  of  this  report  was  to  describe  an  approach  by  which  atomistic  physics  can  be 
embedded  into  a  continuum  formulation  for  large  scale  systems.  This  goal  has  been  achieved 
by  formulating  a  consistent  set  of  equations  involving  a  classical  atomistic  potential  at  the 
fine  scale  and  general  finite  strain  and  deformation  elasticity  at  the  coarse  scale.  Simple 
I'D  analytical  results  and  2-D  numerical  experiments  were  shown  to  illustrate  the  approach 
and  its  features.  More  realistic  multiaxial  problems  in  two  and  three  dimensions  for  more 
detailed  validation  are  the  subjects  of  ongoing  work. 
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